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ment under Markovian approximation, and a derivation of the uncertainty principle at finite 
temperature for a composite object, modeled by two interacting harmonic oscillators. The 
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entanglement, fluctuations and dissipation of mesoscopic objects towards the construction 
of a theoretical framework for macroscopic quantum phenomena. 
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I. INTRODUCTION 

Macroscopic quantum coherence phenomena (MQP) manifested in double slit experiments, 
micromechanical resonators, Bose-Einstein condensates, Josephson junction circuits, mesoscopic 



systems, or even mirrors (see, e.g 



■Hjinstem condensates, Josepnson junction circuit 

.aaaaaaaaa no, mi, im 



Il5l |) is a subject 



of both basic theoretical and practical application interest. Theoretically it focuses on issues at 
the intersection of two trunk lines of important inquires in physics: the relation between the 
microscopic and the macroscopic world on the one hand, and the relation between the quantum 
and the classical on the other. Rapid recent advances in precision measurements with high degree 
of control and adaptability in atomic-optical, electro-mechanical, opto-mechanical, nano-material, 
magnetic-spin and low temperature systems have provided the rationale and substance for such 
theoretical investigations, and in some emergent areas where high goals are set, such as the quest 
for quantum information processing, even with some sense of urgency. 

The issues of interest in MQP include quantum dissipation, entanglement, teleportation, de- 
coherence, noise, correlation and fluctuations. A familiar model which one could use to address 



many of these issues is the quantum Brownian motion (QBM) and its dynamics 



described by the master equation or the associated Langevin or Fokker-Planck equations. But since 
the systems of interest to MQP necessarily involve many microscopic or mesoscopic constituents, 
a many-body generalization of QBM is needed. In addition, since most of these systems involve 
non-negligible correlations amongst their components, quantum memory (non-Markovian) effects 
cannot be ignored. Even for the well-studied single harmonic oscillator (1HO) QBM, Markovian 
approximation is valid only for a high temperature Ohmic bath [17j |. Fortunately an exact mas- 
ter (HPZ) equation [lij for the 1HO with bilinear coupling to a gene ral environment has been 



191 ] and Wigner function [2 



found via several techni ques ranging from the influence functional 
to quantum trajectories [21[. The 1HO master equation for the QBM is complex enough to en- 
compass non-Markovian dynamics yet simple enough to yield exact solutions. (See, e.g., [22] and 
references therein.) The new challenge is to find the master equation for iV oscillators in a general 
environment good for the analysis of these issues in mesoscopic physics. 

In this paper we show the derivation of such an equation for two coupled harmonic oscillators 
(2HO). A key observation is that this problem can be mapped into that of a single harmonic 
oscillator in a general environment plus a free harmonic oscillator. Since the master equation with 
all its coefficients for the 1HO QBM is known [2^] one can derive the master equation for the 
2HO QBM easily from them. As an application of this model, we can deduce the decoherence 
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properties of the 2HO system following the similar pattern of the 1HO. As another example, we 
show explicitly how, in some parameter choice, under the Markovian limit, an entangled state 
evolves into a separable state in a finite time. 

The results derived in this paper may be deduced by intuitive reasoning, but we are not aware 
of any theoretical study which yields our results. Our aim here is to provide a proof, or at least a 
plausibility argument, to the effect that the center of mass coordinate is the one most sensitive to the 
environmental influence. This model and its generalization to N harmonic oscillators are expected 
to be useful for the analysis of quantum coherence, entanglement, fluctuations and dissipation of 
mesoscopic and macroscopic objects. 

The paper is organized as follows: in Section [TH we consider the dynamics of two harmonic os- 
cillators coupled to a common heat bath. By employing the center of mass and relative coordinates 
we show how to derive the master equations of the two coupled Brownian particles. In Section ITO1 
we use the influence functional method and derive an exact form of the propagators for the reduced 
density matrices. These results are expected to be useful for analyzing general statistical mechan- 
ical properties of quantum open systems. In Section [IV] we give three examples as applications of 
this master equation: the quantum decoherence and disentanglement of two interacting Brownian 
oscillators in a general environment, and the uncertainty relation at finite temperature for a com- 
posite object modeled by two interacting oscillators. In Section|V]we mention a few more problems 
and physical issues where the results from this work can be usefully applied to for their analysis 
and further extension of the present study. Technical details are relegated to the two appendices. 



II. THE MODEL AND THE EXACT MASTER EQUATION 

Quantum Brownian motion (QBM) of a damped harmonic oscillator bilinearly coupled to a bath 
of harmonic oscillators has been studied for decades, notably by Feynman- Vernon and Caldera- 
Leggett using path integral techniques [l7]. For such a model an exact master equation can 



be deduced without making the Markovian approximation [191 ] . The purpose of this section is to 
extend the well-known Brownian motion model into the case where the system of interest contains 
two coupled harmonic oscillators. 
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A. The Model 

The Hamiltonian of the total system consisting of a system (sys) of two mutually coupled 
harmonic oscillators of equal mass M and frequency f2 interacting with a bath (bath) of Nb 
harmonic oscillators of masses m n and frequencies tu n in an equilibrium state at a finite temperature 
T can be formally written as, 

Htot = H sys + -ffbath + Hint, (1) 

where 

P 2 1 P 2 1 

H svs = — i- + -Mft 2 x? + —2- + -Ml^x 2 , + k(x! - x 2 ) k (2) 
y 2M 2 1 2M 2 2 v 

is the system Hamiltonian for the two system oscillators of interest, with (xi,x 2 ) displacements, 
conjugate momenta (Pi,F^) and coupling constant re, 

N b 2 i 

#bath = + o 771 "^ 9 ^ ^ 

n=l n 

is the bath Hamiltonian with displacement g n for the n th oscillator and conjugate momentum p n 
and 

N B 

Him = {x\ + x 2 ) ^2 C nQn (4) 
n=l 

is the interaction Hamiltonian between the system and the bath. Here for simplicity, we have 
assumed that the two harmonic oscillators are coupled with the same coupling constants C n to the 
bath oscillators. 

Our primary focus in this paper is to derive an exact master equation for the two coupled 
harmonic oscillators. Since the two harmonic oscillators interact with a common thermal bath, 
there will be induced coupling between the two harmonic oscillators even when initially they are 
uncoupled. Thus, the master equation for 2HO QBM is not simply the addition of the two master 
equations for 1HO QBM. It must account for the mutual interactions between the two Brownian 
particles introduced by their coupling to the common heat bath. Of interest is a comparison with 
the model that consists of 2HO each in its own heat bath. In our model, the coupling to a common 
heat bath can give rise to several new features, of particular interest here is the generation of 
entanglement between the two Brownian particles due to the back-action of the heat bath on the 
system 
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However, as is well-known for classical mechanics, the dynamics of an N body quantum open 
system can be made simpler by changing the N body coordinates to that of their center of mass 
(cm) and relative (rel) coordinates. Here, the difference is that the N harmonic oscillators (NHO) 
are coupled with an environment and we seek a quantum mechanical treatment. A quantum 
mechanical theory of N body dynamics forms the theoretical basis for treating MQP. In this paper 
we treat the 2HO case. We will show in what follows that the exact master equation for the 
two coupled harmonic oscillators can be obtained directly from the master equation for the single 
harmonic oscillator, known as the Hu-Paz-Zhang (HPZ) master equation. 

Let us first rewrite the total Hamiltonian in terms of a set of new variables X,x,P,p defined as 

X = ^(xi + x 2 ), x = xi-x 2 , (5) 

P = P 1 +P 2} p= l -{P 1 -P 2 ), (6) 

and the new masses M\ = 2M, M 2 = M/2. In terms of these new variables the Hamiltonian ([1]) 
takes the following form: 

H sys = H cm + H IC \ (7) 



where 



p2 i 

H cm = + -Mi^ 2 X 2 , (8) 

cm 2Mi 2 V ' 



ffrel = + \M 2 U?X 2 + KX\ (9) 



and 



N B N B N b 

H int = {x 1 + x 2 )^C n q n = 2X^C n q n = X^Cnqn (10) 

n=l n=l n=l 

where C n = 2C n are modified coupling constants. Since ([5]) and ([6]) are canonical transformations, 
all the commutators are preserved, and it is easy to check that 

[X, P) = [x, p) = ih, [P, x] = [p, X] = [X, x] = [P, p] = 0. (11) 

We see that the fictitious particle with mass M 2 and dynamical variables x,p has no interaction 
with either the cm particle with mass Mi with canonical variables X, P or the oscillators of the 
heat bath with canonical variables q n . 
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The total Hamiltonian H to t in (1) can now be written as -f/tot = H[ ot + H re \ with a new effective 



total Hamiltonian 



tot — ^cm + Hmt + -ffbath 
p2 1 



2Mi 2 



N B N B 2 



(12) 



n=l n=l 

This Hamiltonian is formally the same as the Hamiltonian for the single harmonic oscillator in cm 
variables (X, P) coupled to the heat bath with coupling constants C n . Note that for this case the 
spectral density I {to) is given by: 



iVfl g 2 
I(uj) = IT o _" ■ (5(CJ - U n ), 



n=l 



2m„uj r , 



(13) 



which differs from the original spectral density I(uS) by a numerical factor 4. 



B. Density Matrix 



We now consider the dynamics of two coupled harmonic oscillators interacting with a common 
heat bath. The density matrix p evolves in time under the unitary operator: 



p(t) = exp 



.-fftot* 



p(0) exp 



. H tot t 



n 



(14) 



From (|T2|) . it is easy to see that this evolution can be decomposed into two parts, a dissipative 
evolution of the center of mass system, 



p(t) = exp 



' h 



p(0) exp 



' h 



and the unitary evolution of the free harmonic oscillator with mass M\ , 



p(t) = exp 



. H TC \t 

n 



p(t) exp 



. H rc \t 



(15) 



(16) 



where H re \ is the Hamiltonian for the 1HO system with reduced mass M2 = M/2 and x,p variables: 

- 2 1 



H, 



rel 



+ -M 2 n 2 x 2 + Kx k . 



(17) 



2M 2 2 

For technical simplicity we make the usual assumption that the initial state of the total system 
is uncorrelated, 



p(0) = Psys(O) X /9bath(0), 



(18) 



and that the heat bath is in a thermal equilibrium state at temperature T. 
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C. Exact Master Equation 



If we are interested in the detailed dynamics of the system but only the coarse-grained effect of 
the bath we can work with the reduced density matrix obtained by tracing p, the density matrix 
of the total system described by (H|) , over the bath variables [28|, [29J : 



Pr = Tr b athP(i)- 



The reduced density operator for the center of mass system is obtained in a similar way, 



(19) 



Pr = Tr bath p(t). 



(20) 



where p defined in (|15p is the density operator for the effective total system (|12p . The relationship 
between p r and p r is given by 



p r (t) = exp 



. H IC \t 



n 



p r (t) exp 



. H xe \t 



(21) 



Tracing over the heat bath variables in (|15[) leads us to a HPZ type master equation for the 
center of mass variables X, P: 

Pr = ^[H cm ,p r ] + ^-IX^Pr] + ^[X , {P, p r }] + ^-[X, [P,~Pr]] - ^[X, [X,f> r ]]. (22) 

Note here that H cm defined in ([9]) is the Hamiltonian for the center of mass variables X, P only. 
This is the exact master equation for X, P interacting with a thermal heat bath with the spectral 
density I(u>) rather than I (to). As a consequence, the coefficients a,b,c,d in the above master 
equation satisfy the same types of equations given by [n| (or [20]]). only the coupling constants 
and mass are different here. 

From the evolution equation (I16p . the required master equation for the reduced density matrix 
p r (t) is thus obtained, 



^[H ByS , Pr ] + a ^l[X\pr] + b MiX,{P,Pr}] + [P,Pr]\ ^ 



[X,[X,p r ]]. (23) 



The only difference between Eq. (|23p and Eq. (|22p is that the unitary evolution is modified by the 
fictitious harmonic oscillator x,p. 

In terms of the original variables x\, x 2 , Pi, Pi, we get 

1 , ... «!/). ■■> - b(t) { 

d(t). 



ifi [H S y S ,Pr) + + X 2) 2 ,Pr] + -^l X l + X 2,{ P 1 +P2,Pr}] 



2h 2 



+ T^\ x l + x 2, [Pi + P2,Pr\] ~ 4^2"[ X 1 + X 2, [ x l + x 2, Pr])- 



(24) 
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This is the exact master equation for the two coupled harmonic oscillators. In the coordinate 
representation, 

p r (xi,X2,y 1 ,y 2 ) = (xi,X2\p r \yi,y 2 ), (25) 
the master equation can be easily written as: 

2M\dx\ dy\ + dx\ dy\ 
1 " .1 



^ = -2M {^r^ + ^~wJ pr+ 2 Mn2{xl - y!+4 - yl)Pr 



+ -M5Q (t)(xi ~yi + X 2 - 2/2 ) ^ i x l +V1+X2 + y2)Pr 



r 



-iWE(t)(xi -yi + x 2 - y 2 ) p. 

+hA(t)( Xl -yi + x 2 - y 2 ) + ^- + ^- + ^-)Pv (26) 

\dx 1 dyi dx 2 dy 2 / 



A set of new notations in (|26|) is introduced to facilitate easy adoption of results from [19]. In 
particular, 

a(t) = M5n 2 (t), b(t) = 2T(t), (27) 

c(t) = A(t), d(t) = E(t). (28) 

It is often useful to use the Wigner function defined in phase space, which is related to the 
reduced density matrix p r in the following way: 



I U l U 2 U\ U 2 \ 

x Pr [xi - —,x 2 - —;xi + —,x 2 + —,t) . (29) 

In correspondence with (|26|) the Wigner function satisfies a Fokker-Planck equation: 

dW / Pi dW , - dw\ 

dt £<\Mdxi dp) 

+ + „) (A + A) * + 2r(i) (A + jL) [(Pl + p 2) ^ 

9 9 \ 2 - . , . / d d \ ( d d 



+ E(f) [m + m) w + A <" {m + m) \m + m) w - < 30 » 

The time-dependent functions 5£l 2 (t), T(t), A(i), E(t) are derived following the same method used 
by HPZ which can be found in Appendix A. 5. 

In deriving the exact master equation we assumed that the initial state for the two harmonic 
oscillators is a product of a function of the relative coordinates and a function of the center of mass 
coordinates. However, it can be easily shown that the derivation is valid for an arbitrary initial 
state of the system regardless of the condition of separability. 
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D. Markov Approximations 

The derived master equation (I26p is exact, so it is valid in both the Markovian and the non- 
Markovian regimes. Memory effects due to the environment is encoded in the time-dependent 
coefficients. In the high temperature ohmic bath limit, the coefficients become constants and the 
spectral density has the form: 

I(u)=M 1 'yuexp(~Y (31) 

where A is a cut-off frequency. In the so-called Fokker-Planck limit (ksT S> hA), we have 

"00 = 2Mlk h BTl S(s), V (s) = M 7 ^( S ). (32) 

Hence, 5Q 2 = — 2^/5(0), T = 7, A = 0, T, = 2M17/C5T '. The constant coefficients obtained for such 
a model give rise to a Markovian master equation. The Wigner function for the center of mass 



coordinates obeys the Fokker-Planck-Markov equation 



dw cm _ p dWcm Mi ^aw t 



dt Mi dX dP 

d 2 

+ 2M ll k B T—W cm , (33) 

where M x = 2M and ft' 2 = ft 2 + 5Q 2 . 

III. THE INFLUENCE FUNCTIONAL METHOD 

In the last section we showed a simple derivation of the master equation for the reduced density 
matrix and the Fokker-Planck equation for the Wigner function. In general it is difficult to get 
a general analytical solution of the master equation. It turns out that in some cases of interest, 
one can get analytic solutions of the master equation through the influence functional method 



311 ] . Using this method, we can get the evolution operator for the reduced density matrix or the 
evolution kernel for the exact master equation which will be very useful for the study of quantum 
decoherence and disentanglement problems. 

Because of this, in this subsection, we will outline the key steps in the derivation of the master 
equation ([26]) via the path integral method. 

As before, the density matrix of the total system at any time t can be written as 

p(t) = e~ i3 ^ pioy^ir 1 . (34) 
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The reduced density matrix of the system is evolved by the propagator J r from time t = to t as 
p r (x 1 ,x 2 ;yi,y2,t) = [ dq n (x 1 ,x 2 ,q n \p(t)\yi,y2,qn) 

dx dy J r (xi,X2, J/i,2/2,*; x w , x 2 o, j/io, J/20, 0) 



xp r (xio,x 2 o; J/10, J/20; * = 0), (35) 

where we have used the collective notation dx^dyo = dxiodx 2 odj/io<%20- 

The evolution propagator J r can be written in a path-integral representation as 

J r (x 1 ,x 2 ,yi,y2,t;x' 1 ,x' 2 ,y' l ,y 2 ,0) 

2 rxkf rvkf i i 

= 11/ Vx k / ^J/fcexp(-5 5 [xi,x 2 ] - t5 , 5[j/i,J/2]) x J r [x 1 ,x 2 ,j/ 1 ,y 2 ], (36) 

where !F{x\,X2, J/i,J/ 2 ] is the Feynman- Vernon influence functional defined by 

dq' n dq' n dq n p bath (q' n ,q' n ,0) Vq n I Vq n exp{-(S I [x 1 ,x 2 ,q n ] - 
Si[yi,y2,q n ] + S B [q n ) ~ S B [q n })} 
= exp{-(S IF [xi,x 2 , j/i, j/ 2 ])}, (37) 

where S/i? is the influence action. For the QBM model we are considering here, the influence action 
can be written as: 

Si F [x 1 ,x 2 ,yi,y2] = -2 / dsi / ds 2 [Ai(si) + A 2 (si)]r/(si - s 2 )[Ei(s 2 ) + S 2 (s 2 )] 

io io 

+ i / d Sl [ 1 ds 2 [Ai(si) + A 2 (si)]i/(si - s 2 )[Ai(s 2 ) + A 2 (s 2 )], (38) 

where 

Si = ^(xi + j/i), S 2 = i(x 2 + j/ 2 ), Ai=xi-j/i, A 2 = x 2 -j/ 2 . (39) 

Note that the integrand in Eq. (|36|) is Gaussian, hence the integral can be computed exactly and 
the explicit form of J r is, 

J r = NexpC-S!-S R ), (40) 

where the expressions of Si and Sr can be written in more compact forms with the following 
notations: 

xt = x\k + x 2 k, yt = yik + V2k, ( 41 ) 

X^ = Xi k -X 2 k, J/fc = J/ifc - J/2fe, (42) 
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and 



whence 

Sj = b x {x+ + yt)(xj - yt) + h(x+ + y£)(xt - yf) 

- h{xf + yt)( x o ~ Vo) ~ b *i x o + Vo)i x o ~ Vo) 
+ h(xt + Vt)(xt ~Vt) + h &( x o + Vo)( x t ~ Vt) 

- b 7 (xt + yt)(xQ - 2/0") - 6 8 (^ + Vo)( x o ~ Vo)> ( 43 ) 

Sr = au(xf - yf) 2 + a 22 {xQ - y^f 

+a 12 {x+ - y+)(xf - yt). (44) 

The functions bi(t) and Oy(i) depend on the environment and can be constructed from the solutions 
to the equations 

b 2 (t) = hi 1 {t), h(t) = i« 2 (f), 6 6 (t) = ^!(t), b 5 (t) = ^w 2 (t), 
& 4 (t) = 1^(0), b 3 (t) = iu 2 (0), 6g(t) = ^(O), b 7 (t) = ^w 2 (0), (45) 

where Wi(t) are functions which satisfy the following equation 

E(s) + 2 S(s) = 0, (46) 

with the boundary conditions: 

Wl (0) = l = w 2 (t), Wl (t) = = w 2 (0), (47) 

a ij( t ) = 7;[ ds i [ ds 2 Ui(si)v{si - s 2 )uj(s 2 ). (48) 
* Jo Jo 

With the expression of J r , we can derive the master equation for the reduced density matrix 
(1261) , This is shown in Appendix A. 

An exact form of the evolutionary operator for the reduced density matrix is a priced object: 
Not only can one derive from it the exact master equation for the reduced density matrix, with 
this explicit expression of the evolutionary operator, given any initial reduced density matrix p r 
at time to one can calculate p r at any later time t without having to solve the complicated second 
order partial differential equation with time-dependent coefficient functions. 

For example, we will apply this evolutionary operator to the study of the decoherence and 
disentanglement of two coupled harmonic oscillators in a common heat bath. One can also use 
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it to calculate the higher moments of physical observables of interest such as the position and 
the momentum operators which enter into the derivation of a generalized uncertainty principle 
for composite objects at finite temperature .32]. It can also be used to address the issue of the 



influence of entanglement on the relation between the statistical entropy of an open quantum system 
and the heat exchanged with a low temperature environment such as studied in [jjjj]. Another 
interesting application would be the entanglement between a qubit and an oscillator. Adopting 
a level reduction scheme, Shiokawa and Hu 



411 ] used the evolutionary operator of 1HO QBM to 
study the dynamics of the spin-boson model. The explicit expression of the evolutionary operator 
for the 2HO QBM may be used to construct effective lHO-spin-boson models found in many 
condensed matter quantum computer schemes for the analysis of the interaction between a qubit 
and a harmonic oscillator and their decoherence and disentanglement dynamics in the presence of 
a general environment. See Section |V] for a more detailed exposition of further applications and 
extensions. 



IV. APPLICATIONS: QUANTUM DECOHERENCE AND DISENTANGLEMENT, 
UNCERTAINTY RELATION FOR A COMPOSITE OBJECT 

In this section we give three examples for the application of this master equation: the deco- 
herence and disentanglement of two coupled harmonic oscillators in a common heat bath, and a 
derivation of the uncertainty relation at finite temperature for a composite object modeled by two 
harmonic oscillators in a general environment. For some simplified cases we obtain analytic results 
which show interesting features such as finite-time disentanglement 



A. Dynamics of Quantum Coherence 

We will assume that the system and the environment are initially uncorrelated. The total 
density matrix at time t = then factorizes into a product of density matrices for the system 
and the environment. As usual, we further assume that the environment is initially in thermal 
equilibrium at a given temperature T. 

We assume initially the 2HO (labeled as 1 and 2) are separated with distance 2Lq and the initial 
wave function of the 1-2 system is given by 



*(xi,x 2 ,t = 0) = s 1 *i(s 1 )*i(s a ) + S2*i(a;i)*2(a;2) 

+s 3 ^ 2 (xi)^i(x 2 ) + 54*2(^1)^2(3:2), 



(49) 
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where we have defined the displaced Gaussian states as 

*i,2 0*0 = N exp [- (X ^° )2 ] exp (±iP Q x), (50) 

and Si are any complex numbers subject to normalization conditions. (We use 1,2 to label different 
initial positions of the center of the Gaussian wave function of harmonic oscillators while x, y label 
different time paths.) 

With an initial reduced density matrix 

Pr(x w ,x 20 ;y 10 ,y 20 ;t = 0) = (x w ,x 20 | *(0))(*(0) | 3/10,3/20) 

= ^2 s^* pij (x 10 , x 2 o; yw, y2o; t = 0), (51) 

the reduced density matrix at t is given by 

p r (x 1 ,x 2 ;y 1 ,y 2 ;t) = J dx Q dy Q J r (x 1 ,x 2 ,y 1 ,y 2 ,t;x 10 ,x 20 ,y 10 ,y 2Q ,0) 

xp r (x w ,x 20 ;y 10 ,y 20 ;t = 0). (52) 

Because the QBM model is linear and the initial state is Gaussian, we can solve the master 
equation exactly for the dynamics of the 2HO system interacting with an environment with a 
general spectral density at any temperature. Therefore, we can obtain the total density matrix if 
the explicit solutions for each component are known, 

Pij(xi,x 2 ;yi,y 2 ;t) = J dxody J r (x 1 ,x 2 ,y 1 ,y 2 ,t;x 1 o,x 2 o,y w ,y 2 o,0) 

xpij(x w ,x 20 ;y 10 ,y 20 ;t = 0). (53) 

Note that since J r and pij are in the form of an exponential with an exponent which is a 
quadratic function in (xio, x 2 o; 2/10j 3/2o)> we can use a standard trick for the evaluation of the 
integral, 

Pij(t) = J dx dy Jt x pij(t = 0) 

= J dx dy exp \-x T ■ Gij ■ x + ^F? ■ x + ^x T ■ + aj] 

C./^ 4 1 -1 — » 

^\nr,j I ■ ('-,f ■ /•;,)• (54) 



y/det 

where x T = (x w , x 20 , 3/10, 3/20) ■ 

Once we have pij(x\,x 2 ; 3/1, y 2 ; t) we can perform the following substitution x\ ^> X\ — 4f; x 2 1— > 
X 2 — tt ; 3/1 i—*- X\ + ; 3/2 X 2 + 4?- and then do the Fourier transform to get the Wigner function 
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at a later time t: 



W l! (X l .X 2 . r\.P 2 .l) = I I ^0 exp (iP lZl + iP 2 z 2 ) 



xptjiXi f^i + + f ;*)• (55) 

Since after the substitution the exponent of pu is quadratic in Zi,Z2, the above integration can 
be evaluated explicitly. These solutions (|54p and (|55p will be useful in decoherence and disentan- 
glement analysis below. The detailed results and the explicit expressions of pij can be found in 
Appendix B. 

When viewed from the center of mass coordinate the physics of decoherence for a 2HO system 
is essentially similar to that described in 



for 1HO 

because the environment couples to the system only through the center of mass coordinate X and 
is independent of the relative coordinate x. The evolution of the relative coordinate part in the 
reduced density matrix is unitary and hence will not affect the decoherence processes. One can 
easily recognize these features from (122p and (|2ip . The effects of environment-induced decoherence 



are encoded in the coefficient functions a(t), b(t) , c(t) , d(t) of (|22p . As one can see from this example 



35] 



four of the matrix elements pn, pi4, pa, Pm are similar to those in the example considered in 
sans the relative coordinates. 

However, the issue of disentanglement is quite different because usually the entanglement mea- 
sure is related to the global property of the whole reduced density matrix. In general, entanglement 
involves both the center of mass and the relative coordinate dynamics. It is difficult to make any 
prediction on how disentanglement evolves from the information of only the 1HO system. For 
instance, while the cm coherence always disappear asymptotically, in contrast, entanglement of the 
two particles may terminate in a finite time. In the third subsection, we will address this issue 
with a simple illustrative example. 



B. Uncertainty Principle for Composite Objects 

In this subsection, the generalized uncertainty relation for a composite object is investigated 
from the viewpoint of quantum open systems. Here the system is modeled by two harmonic 
oscillators and the environment by a heat bath at temperature T. As such, both thermal fluctuation 
and quantum noise come to play when the uncertainty relation between position and momentum 
is considered 
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The exact solution for the two harmonic oscillators coupled to a common heat bath can be found 



15 



by decomposing the total system into two fictitious surrogate subsystems, namely, the subsystems 
described by the center of mass and the relative coordinates, respectively. Such a decomposition 
guarantees that the two subsystems are decoupled, and as such, the solution of the total system is 
a tensor product of the two subsystems: 



Pern ® /Orel- 



(56) 



Using the center of mass coordinate as described by the Hamiltonian (pQ) , the complete information 
about the state of the open system is contained in the reduced density operator p r (t). 
For a class of initial Gaussian states given by 



ijj(x, 0) = iVo exp 



(x - xpf 
4a 2 



+ tPox 
n 



(57) 



where, Nq = l/(27r<7 2 )4, the initial density operator for each fictitious harmonic oscillator in the 
coordinate representation can be written as: 

(x — Xq) 2 (x' — Xq) 2 



p(x,x',0) = ip*(x,0)ip(x',0) = iV 2 exp 



4a 2 



4a 2 



-^Pox + -pox 



(58) 



In order to compute the variance of position and momentum operators, it is more convenient to 
use the Wigner function which can be obtained from the Fourier transform of the density operators 
(|29j) . To be more specific, for the harmonic oscillator representing the center of mass degree of 
freedom, the corresponding Wigner function is simply given by: 



W cm (X,P) =iV 2 exp 



2(J 2 - W ( P -P0) 



(59) 



The variance of the operator X,P denoted by (AX) 2 = (X 2 ) - (X) 2 and (AP) 2 = (P 2 ) - (P) 2 
can be computed easily 

1 



(X 2 ) 
(P 2 ) 



2irh 
1 

2ttE 



dXdPX 2 W cm (X,P,t) 



l — j dXdPP 2 W cm (X,P,t) 



(60) 
(61) 



where W cm (X, P, t) is the solution of the Fokker-Planck equation for a single harmonic oscillator 
(see Appendix B or 32j|). In particular, for an ohmic environment (|32p . the uncertainly relation in 
the weak damping limit (7 << 20) is given by 



U(t) = (AX) 2 (AP) 2 (Ax) 2 (Ap) 2 = f cm (t)f iel (t) 



(62) 
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with 



h 2 

/cm (t) = ~T 



e -7t + coth _^_ (1 _ e -7t) 



+ TTcoth 



Ml' 



/rel(t) 



+ ft 2 

4 



2k B T 
1 - 5 2 



2k B T 

(l-^) 2 
45 



;i- e -^)- (1 ^! )7 siu2(// 



45 



sin + ^ I coth 



20' 



8WS 

1 + S 2 
25 



-it 



sin 2 fl't 



-i-ft 



1 



l + -^r(l-<5 2 ) 2 sin 2 20t 
4<5 2 v ; 



where S)' = y^ 2 - 7 2 /4 and 5 = 2VLo 2 /h. At short times (i « 1/7, 1/0), 



fcrait) — . 



1 + 2 (5 coth 



2/c B r 



1)7* 



4 ' 



(63) 
(64) 

(65) 
(66) 



In this short time span, the time-dependent quantum dispersion of the wave packet constructed in 
the relative coordinate may be ignored. It is interesting to compare the uncertainty relation (|62p 
with that between the x±,x 2 and pi,p2, denoted by U XiPi : 



1 



C/ XiPi = (Ax 1 ) 2 (A Pl ) 2 (Ax 2 ) 2 (A P2 ) 2 > -U(t) 



(67) 



As will be shown in the next subsection, the variance of the operators x and X etc can indeed 
provide some useful information about the evolution of quantum entanglement of the Gaussian 
states. 



C. Dynamics of Entanglement: An Example 

As shown in Subsection IIV Al the decoherent effects of a thermal heat bath is captured by the 
influential functional appearing in (|52p . An environment that destroys quantum coherence can also 
disentangle two quantum Brownian particles. The dynamics of decoherence and entanglement of 
two harmonic oscillators interacting with a common environment is useful for understanding some 
basic issues in macroscopic quantum phenomena. We will present a more detailed study of this 
issue in a later paper. Here we show a simple example which has analytic solutions. Take as initial 
state the Wigner function: 

W {x ! , x 2 , P x , P 2 ) = W cm ( X, P) W rel (x, p) = e~ & ~ ^ e~& " . (68) 

where P, X, x and p are canonical variables defined in (|5|) and ([6]) . We have omitted an irrelevant 
normalization factor. Note that the widths a 2 , b 2 ,c 2 and d 2 cannot be chosen arbitrarily since they 
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have to satisfy the uncertainty relations: 

fr 2 h 2 

a 2 b 2 > 1. c 2 d 2 > IL ( 69 ) 
- 4 - 4 

For a wide range of parameters a, b, c and d, the Wigner function W(X, P, x,p) is entangled, since 
generally it cannot be written as a product of W\(x\, P\) and W9.(x 2 , Pa). At any time £, it is 



known that the separability of the state f)68[) can be easily detected |36l. l37|. 

Now we consider the dynamics of this state under the influence of a common environment. For 
greatest simplicity, we assume two free particles coupled to a Markovian thermal bath (Setting 
£1 = and k = 0) and assume the dissipation in cm coordinates is negligible. Under these 
conditions, the Wigner equation W cm (X, P) for cm coordinates ([33]) takes on a simple form: 

dw cm p dw cm , d 2 w cm 



+ D ^f, (70) 



dt Mi dX dP 2 

where D = 2Mi^ksT. The solution for the dissipative evolution of the center of mass can be easily 
obtained, and from it, we can compute the variances of X and P at time t to be: 

9w s 2Dt 3 b 2 t 2 , , N 

iAX) ^ = up + uP + "' (71) 

(AP 2 )(t) = 2Dt + b 2 . (72) 

Since the evolution of the Wigner function W Te i(x,p) for the relative coordinates x,p is unitary, 

dW rel p dW Te i 



dt M 2 dx 

the variances at t are simply given by 



(73) 



(A* 2 )(£) = ^ 2 + c\ (74) 

(Ap 2 )(t) = d 2 . (75) 

According to [36|], we may choose the EPR-like operators as : 

u = x\ - x 2 , v = Pi + P 2 , (76) 

where Xi,Pi (i = 1,2) are the dimensionless variables satisfying [xj, Pj] = iSij, 

* = (1?) % - # ' = (iii7o) ,pi ' <i = 1 ' 2) - (77) 

Then the Gaussian state (|68p at £ is disentangled if and only if the following inequality is satisfied 

(Au 2 )(t) + (Av 2 )(t) > 2. (78) 
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Inserting ([72"]) and ([74"]) into the above inequality, one gets, 

,4i 2 + St + C > 2, (79) 

where 



4d 2 I MD 



A = iPVur' (80) 



B = 2^. (81) 
b 2 n I MD 



From ([79]) . the disentanglement time tdent can be determined to be 



_ _ B + ^ B 2 _ AAC + 8A 

Wt - 2^ ' ^ ' 



Thus after i > ident the state ([68]) becomes completely separable. 

In situations when the 2HO are coupled or share the same environment, it is expected that for 

ing 



38] 



some initial states entanglement will persist longer than the case when there is no direct coup 
between the two oscillators and each of them is coupled to a separate environment (See, e.g. 
for two qubits in a common electromagnetic field) . This is what one might anticipate would happen 
for our model in the more general cases. On the other hand, as shown in this simplified example, 
finite-time disentanglement may yet occur for some initial states when there is no direct coupling 
between the two oscillators. 

Such finite-time decay behavior has been noted before in several cases where two qubits 
or two harmonic oscillators 



34, 



33] 



391 ] are individually coupled to their own heat baths. We show 
here the onset of the finite-time decay for the case of a common heat bath. However, it should be 
emphasized again that the finite-time disentanglement process found here depends crucially on the 
choice of initial states because for some initial states the mutual actions between the two harmonic 
oscillators may lead to entanglement generation. As shown in the case of two-qubits under phase 
noises, when the initial states are protected by a decoherence-free subspace quantum entanglement 
is shown to be robust against the thermal noise [^(J. The 2HO model considered here will exhibit 
similar features, but further details will go beyond the scope of this paper. 



V. FURTHER APPLICATIONS AND DEVELOPMENTS 



a. Summary In this work we derive an exact master equation for two coupled quantum 
harmonic oscillators interacting via bilinear coupling with a common environment at arbitrary 
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temperature made up of many harmonic oscillators with a general spectral density function. We 
first show a simple derivation based on the observation that the two harmonic oscillator model 
can be effectively mapped into that of a single harmonic oscillator in a general environment plus a 
free harmonic oscillator. Since the exact one harmonic oscillator master equation is available [3] 
the exact master equation with all its coefficients for this two harmonic oscillator model can be 
easily deduced from the known results of the single harmonic oscillator case. In the second part we 
give an influence functional treatment of this model and provide explicit expressions for the evolu- 
tionary operator of the reduced density matrix which are useful for the study of decoherence and 
disentanglement issues. We show three applications of this master equation: on the decoherence 
and disentanglement of two harmonic oscillators due to their interaction with a common environ- 
ment and a derivation of the uncertainty principle at finite temperature for a composite object, 
modeled by two interacting harmonic oscillators. For the example of entanglement dynamics under 
Markovian approximation we find finite-time disentanglement taking place for a Gaussian state. 

b. Decoherence and Disentanglement We mention some further developments and applica- 
tions where our analysis of the 2HO QBM model can be usefully extended to or compared with. 
First, for the study of decoherence and disentanglement between two observers, a direct comparison 
can be carried out with some recent findings in 43] where the model of two harmonic oscillators 



in relativistic motion (one could be in uniform acceleration) in a common field in Minkowsky or a 
black hole spacetime. In the latter situation it is of interest to see how entanglement and teleporta- 
tion will be affected by its unusual causal properties. The case of two oscillators in inertial motion 
in ordinary Minkowsky spacetime would correspond to our problem here after invoking Lorentz 
invariance. Second, pursuant to our analysis of the uncertainty principle for composite objects, the 
substance of our calculations there could be applied to another interesting physical issue pertaining 
to the Landauer principle 



441 ] and the Clausius inequality. Landauer principle which rests at the 
foundation of the thermodynamics of information processing, states that (paraphrased in the words 
of Bennett [3]) "any logically irreversible manipulation of information, such as the erasure of a 
bit or the merging of two computation paths, must be accompanied by a corresponding entropy 
increase in non-information bearing degrees of freedom of the information processing apparatus or 
its environment. Conversely, it is generally accepted that any logically reversible transformation 
of information can in principle be accomplished by an appropriate physical mechanism operating 



in a thermodynamically reversible fashion". (See also 461. 1471. l48j] . the last contains a proposal for 
a generalized Landauer's principle.) It is well known that the root of this relation is the second 
law of thermodynamics, but how to measure a logical operation in physical terms or to associate a 
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logical state or its transformation with an energy cost and an entropy increase is a new challenge. 

c. Quantum Information and Thermodynamics There are many angles to see how Landauer's 
bound in quantum information theory is related to Clausius' inequality in classical thermodynam- 
ics. One such approach is by way of quantum open systems which can treat the dynamics of the 
system and its quantum information content in fully nonequilibrium settings. This is the concep- 
tual framework and technical systematics we have adopted. Here, dissipation and decoherence in 
the system and disentanglement between the system and its environment may be followed closely 
by the evolution of the reduced density matrix (RDM), and the entropy change of the system in 
the thermodynamic limit may be calculated, with little difficulty In this vein, using the quantum 



Brownian model (QBM) of the Caldeira-Leggett (CL) type Hoerhammer and Buettner 



501 ] inves- 



tigated the influence of entanglement on the relation between the statistical entropy of an open 
quantum system and the heat exchanged with a low temperature environment. (See also 491] ). 
Their two Brownian oscillator model is of particular relevance to our work here. Compared to 
the case of a single Brownian particle, two coupled harmonic oscillators can account for how the 
internal degrees of freedom of the system would affect the heat and entropy changes. Because they 
adopted the CL treatment their results are subcases of ours here (in the same way that the CL 
treatment \l\ of QBM is related to the HPZ treatment [lfl ]. viz, the latter preserves the positive 
definiteness of RDM in its entire evolution and the HPZ master equation extends the range of valid- 
ity to non-Markovian regimes.) The CL results are valid only for ohmic baths at high temperatures 
pertaining to the Markovian regime. For low temperatures and nonOhmic baths pertaining to the 
nonMarkovian regimes the HPZ treatment is expected to yield more accurate results. Thus using 
the master equation presented here for the 2HO QBM model following HPZ treatment and the 



analytical solutions found recently 



511 ] for various parameter ranges one could obtain an improved 



Landauer bound for quantum information processing in the nonMarkovian regimes. On the other 
side of the balance, the Clausius inequality, operative only in the thermodynamic limit, would be 
too coarse a measure for the energy cost and entropy change of quantum information processing 
anyway. With the master equations derived here there is much room for tightening the Landauer 
bound. 

d. Qubit - Oscillator Entanglement As subcases of the present model one can investigate the 
interaction between a two-level system with a harmonic oscillator in a general environment which 
is of general interest for quantum computer design purposes. One could apply a level reduction 
scheme such as that used in to one of the two harmonic oscillators, turning the 2HO-bath 
model into an effective lHO-spin-boson model where the bimodal oscillator mocks up a qubit. 
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The simpler case without an environment would correspond to a two level atom in a multi-mode 
cavity, such as studied in 



421 ] . Doing a level reduction scheme for both oscillators and viewing 
the harmonic oscillator bath as a field would reduce our 2HO QBM model to that of two qubits 
interacting either directly or indirectly through a common field. An example of the latter situation 
is studied in [38]. One can use the exact master equations here under appropriate simplifications 
to describe the nonMarkovian dynamics of such systems. 

e. Quantum Superposition of Two Mirrors As mentioned in the Introduction, a new category 
of problems which has received much attention lately is represented by the quantum superposition 
of two mirrors 8]. The two mirrors can be modeled by two quantum harmonic oscillators, but 
in most models for proposed experimental designs, the mirrors are coupled by radiation pressure. 
This class of models with photon number - mirror displacement (Nx) type of coupling used for 
mirror-photon entanglement [521 ] . ent ang lement cooling of a mirror 53] and entanglement of test 



masses and standard quantum limit [54| is very different from the class with bilinear coupling in 



QBM studies (Beware of inconsistencies in the usual master equations for this problem, see [55]|). 
On the surface the convenience of the 2HO model which possesses many useful solutions would 
not be readily available, but a recent observation by Galley could provide a bridge to these two 
common classes of models and unleash the resources gathered from the 2HO QBM problem for the 
solution of this type of quantum optics problems. (See [5a].). 

/. Macroscopic Quantum Phenomena Finally, a whole range of issues in macroscopic quan- 
tum phenomena can be addressed with the master equation (or the associated Langevin or Fokker- 
Planck equations) derived here. In particular, decoherence and disentanglement in 2HO system 



231 ] are currently under study. 



under more general conditions and A-harmonic oscillators systems 
It can also be applied to the analysis of quantum decoherence, entanglement, fluctuations, dissi- 
pation and teleportation of electro- opto-mechanical systems and superposition of moving mirrors 
due to quantum and radiative effects. 
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APPENDIX A: DERIVATION OF EXACT MASTER EQUATION FROM PATH 

INTEGRAL 

Deriving the master equation from the path integral is lengthy, but one of the advantages of 
this derivation is that the explicit form of the propagator can be used to find an explicit solution 
of the equation in many interesting cases. We will mainly follow the steps in [ly] and outline the 
key steps in deriving the master equation from the path integral method. 

From (|35p . it is easy to see that, to get the master equation, one first needs to calculate 
J r {t + dt, 0) — J r {t, 0). The complete derivation can be decomposed into the following four steps. 



1. Step one 

Our first task is to take the functional representation of J r {t + dt, 0) and divide each of the path 
integrals into two parts. We introduce four intermediate points x\ m ,x 2m ,y\ m ,y 2m at time t and 
integrate over them, thus symbolically, we write 

t+dt;xif roo ?t\x\ m rt+dt;xif 

Vx\ I dx lm / Vxi \ Vx x . (Al) 



'0;xio J -oo J0;x\q J t;x\ m 

There are three similar expressions for the sum over X2,y 1,2/2 histories. 

The original histories x\(t) are functions defined on (0, t + dt) time interval with xi(0) = 
xio, x\(t + dt) = xif. The new set of histories xi(r), x\{r) are functions defined on (0, t), (i, t + dt) 
intervals with xi(0) = x w ,xi(t) = (t + dt) = xif. 

So we can write 

A[xi,x 2 ,yi,y2] = S s [xi,x 2 ) - S s [yi,y2\+ 5A[xi,x 2 ,yi,y2\ 

= A[xi,x 2 , yum] + A[xi,x 2 , yum] + A[xu^2,yum,xux 2 ,yum], ( A2 ) 

where Ai term mixes the x histories with x ones. The appearance of the A{ term is due to the 
non-locality of the influence functional. 
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2. Step two 

Next, we will use straight line histories approximation of (xi, x 2 , Vi, y 2 )- First, note that 

Xl(s) = Xi m + (x lf - X lm )?-jj- = X lm + /?l x ^-, ( A3 ) 

and similarly, 

s — t s — t 

X2(s) = X 2m + {x 2 f~ X2m)—jj- = X 2m + 02x~^, (A4) 

s — t s — t 

yi(s)=yi m + Piy —jj- , m{s) =y 2m + fay ■ (A5) 

To compute the time derivative of J r , take the limit dt — > 0. Thus we can write 

2 pt+dt;x k f rt+dt;y k f 



* rt+at;x kf rt+at;y kf ■ 

\[ Vx k ^yfcexp(-^[Xi,X2,yi,2/ 2 ]) 

k=l J°'' x kO J0;y k0 n 

= / daW^fcmexp(-A[£i,X2,yi,y 2 ]) 

JL rt;xkm rt;ykm ^ ^ 

x j Vx k I Py fc exp(-A[xi,X2,yi,y2])exp(-yl4xi,X2,yi,y2,xi,X2,yi,y2])(A6) 

k=1 J0;x k0 J0;y k0 n h 

Expanding A in dt and keeping the contributions of the first order terms, we get, 
A[x u X2,yi,U ~ 2m(Pix + PL ~ Ply - Ply) ~ ^rnn 2 dt(x 2 lf + x 2 2f - y 2 lf - y 2 f ) + ■■-, (A7) 



and 



Ai[xi,X2,yi,y2,xi,x 2 ,yi,y2] ~ -dt / dsJg(s)(Ei(s) + £ 2 (s)) 

Jo 

+idt [ dsJ K {s)(A 1 {s) + A 2 (s)), (A8) 

JO 



where 



o rt+dt 

J Sl + J S2 = J g ( s )_ jf ^'(AxO/) + A 2 ( S ')Ms' - s) 

~ 2(xij — yi/ + x 2 / — y 2 f)r](t — s) H , (A9) 



and 



i rt+dt 

Ja 1 + Ja 2 = <^a( s )^ y ds'(Ai(s') + A 2 (a'))f («' - a) 

~ ~~ 2/1/ + x 2 f — y 2 f)v(t — s) -\ . (A10) 



Here we can keep only terms up to the first order in [3 2 . 
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In summary, the propagator J r can be formally written as 

J r (x lf ,x 2 f,y 1 f, 2/2/,* + dt\x w , x 2 o, ?/io, 2/20, 0) (All) 

/OO f'OO /'OO /'OO 

d/9ix / dft* / ^ / d/3 2 , exp (^r(/?L + #L " ft v ~ Ply)) ( A12 ) 
-oo J — oo J —oo J — oo LuXXl 

x {1 - ^dtp^rci/, z 2/ ) - y(2/i/, 2/ 2 /)]} 

x J r (x lm , x 2m , 2/im, 2/2m, t + di|xi , ^20, 2/io, 2/20, 0; [6]), (A13) 



where 



J r (x lm ,x 2m , 2/im, 2/2m, * + cft|xio, ^20, 2/io, 2/20, 0; [6]) = (A14) 
rt;xi m rt;x2m rt;yi m rt;y2m ^ 

/ Vx\ j Vx 2 / Pf7i / ^2/2 exp (-A[x4,x 2 , 2/1,2/2]) (A15) 

JO;xio ^0;a;20 JO;yio JQ\V20 

exp [l(-di / dsJ g (5)(Ei(a) + E 2 (s)) + idt [ dsJ a (s)(Ai(s) + A 2 (s)))], (A16) 

" J0 JO 



and 



(A17) 



where the sources 6 are functions of the end points. Note that J r (b) can be interpreted as the 
evolution operator under the action of two external sources. 

3. Step three 

Computation of the path integral J r (o) can be done as follows. First, one can show that 

Jr(x lm ,x 2m ,y lm ,y 2m ,t\x w ,x 20 ,y w ,y 20 ,0;[b]) = (A18) 
J r (xim,x 2m ,y lm , y 2m ,t\x w , x 20 , 2/10, 2/20, 0)W(x lm , x 2m ,y lm , y 2m , x 1Q , x 20 , y w , y 20 , di)(Al9) 

(Note that the function J r is the evolution operator without source while the function W is a 
function of the end points. ) 

Then we may parametrize the paths, and write 

Ei(s) = ^i(s) + S cU (s), E 2 (s) =v9 2 (s) + S c/i2 (s) (A20) 
^(^(sJ + AciW, A 2 (s)=^(s) + A cJ>2 (s) (A21) 

where the "classical paths" I \ are the solutions to the equation of motion derived from the 



A 



real part of A[Ei, E 2 , Ai, A 2 ]. 



cl 
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After this path reparametrization and making a saddle point approximation, this path integral 
J r (x lm , x 2m , yim, y 2m , t\x 10 , x 20 , yio, 2/ 2 o, 0; [b]) can be written as 



J r (xi m ,x 2m ,yim, V2m, t\x 10 , x 20 ,yi , 2/ 2 o, 0; [b]) = J r (0, 0, 0, 0, t\0, 0, 0, 0, 0; [b]) 
x exp (-A[T, d;1 , S d)2 , A d;1 , A dj2 ]) 



where 



x exp[-(-dt y dsJ f (s)(S cM ( S ) + S d)2 ( S ))+^ y dsJ A (s)(A djl ( S ) + A dj2 (s)))],(A22) 



rt;ipi=0 ft;<P2=0 rt;xpi=0 rt;ip2=0 

J r (0,0,0,0,t|0,0,0,0,0;[6]) / 2tyi / Vip 2 / P^i / 2^2 

J0;ipi=0 J0;(p2=0 J0;V>1=0 J0;i/> 2 =0 



exp{z[/'ctei r^ T ( S i)6( S i, S2 )^(s 2 )+ f ds^ T (s)-B(s)}}. 
Jo JO 2 Jo 



(A23) 



Note that 



*2 



and 



/ _ 



B 



V?2 



dtJ^ \ 



(A24) 



idtJ^ + i 
—dtJ^ 



(A25) 



where is a new source which appears because the nonlocality of the influence functional. It 
couples the classical paths to the ^ paths. 

rt 



J K (s)= [ ds'[A dA (s') + A cl ^s'Ms-s'). 
Jo 

The matrix operator 0(s\,s 2 ) is defined as follows: 



(A26) 



On(si,s 2 ) = 3 3(-si,S2) = Oi 3 (s 1 ,s 2 )0 31 (si,s 2 ) = 0, 



(A27) 



2 2(si, s 2 ) = O44O1, s 2 ) = 24 (si, s 2 )0 42 (si, s 2 ) = 2zi/(si - s 2 ), 



(A28) 



Oi4(si, s 2 ) = 32 (si,s 2 ) = 26(s 2 - si)r)(si - s 2 ), 



(A29) 
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41 (s u a 2 ) = 23 (si, a 2 ) = 20(«i - s 2 )»/(«i - 82), (A30) 
12 (Sl,S 2 ) = 34 (S1, S 2 ) = (J^~2 + fi 2 ^ 5(8! - 8 2 ) + 20(s 2 - si)7?(si - s 2 ), (A31) 

2 i(si, a 2 ) = 43 (si, s 2 ) = (-^2 + ^1 5 ( s i " s 2) + 20(«i - s 2 )r?(si - s 2 ). (A32) 



The Gaussian path integral can be computed in terms of the inverse of the operator O, which 
is given by G = O^ 1 . Hence to first order in dt, we have 

rt;ipi=0 rt;<P2=0 ft;ipi=0 rt;ip2=0 

J r (0,0,0,0,t|0,0,0,0,0;[6]) / Vy> 2 / Xtyi / % 

J0;<pi=0 Jo;ip 2 =0 J0;ip!=0 J0;^ 2 =0 

exp{z[/*d Sl r^ T ( Sl )0(si, S2 )^( S2 )+ f ds^ T (s)-B(s)}} 

JO JO 2 JO 

= Jv<PiJ V( P2 J Vih j V^ 2 exp {^(# T + £ T • 6 _1 )0(^ + O" 1 • 5) - ^iFO^S]} 
= Zo^expj-^^O- 1 ^} 
« Z (*)(l - l -B T 6- x B) 



k, Z (t)(l-dt / dsi / ds 2 Jg(si)[Gi 2 (si,s 2 ) +Gi 4 (si,s 2 ) + G 2 i(s 2 ,si) + G 4 i(si,s 2 )]J 5 (s 2 )). 

Note that the Green's function (Gi 2 + G 32 ) = Gi 2 (si,s 2 ) satisfies the following equation 
d 2 f si 

^|Gi 2 (si, s 2 ) + ft 2 Gi 2 (si, s 2 )+4j drriis! - t)G 12 { Si ,t) = 5( Sl - s 2 ) (A33) 

with boundary conditions Gi 2 (0, s 2 ) = G\ 2 (si,t) = 0. The equation for (G21 + G 23 ) = G 2 \(s\,s 2 ) 
are analogous. 

Now we can show that 

J r (xi m ,x 2m , yim,y 2m , t\xi ,x 20 , yw,y 2 o, 0; [b]) 

= J r (0, 0, 0, 0, t|0, 0, 0, 0, 0; [6]) exp {i(A[X d>1 , E dj2 , A djl , A d , 2 ])} 

xexp{i(-dt / dsJg(s)(E cJi i(5) + E c , i2 (s)) + idt / dsJ^(s)(A c/jl (s) + A d;2 (s)))} (A34) 
Jo Jo 



Z (t) exp {M[E d)1 , £ d;2 , A di i, A d , 2 ]} (A35) 
'o Jo 



:{l-dt [ ds! [ ds 2 Jg(si)[Gi2(si,s 2 )+G2i(si,s 2 )]JA( s 2) (A36) 
Jo Jo 

-idt [ dsJ^s)(Z dA (s) + Z d , 2 (s)) + (i) 2 dt ( dsJ K (s)(A clA (s) + A cl)2 (s))} (A37) 
Jo Jo 

= Jr(xim, x 2m , yim,y 2m , t\x w , x 20 , y w ,y 2 o, 0) 

xW(x lm , x 2m , y lm ,y 2m ,x w , x 20 , y w ,y 20 ,dt), ( A38) 
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where W is given by, 

W(x lm , X 2m ,yi m , V2m, ^10, ^20, 2/10, 2/20, dt) 

= l-idt[f ds2(A lf + A 2f )v{t - s) Ul (s)t d (0) 
Jo 

+ [ ds2(A lf + A 2f ) V (t - s)u 2 (s)t cl (t)} (A39) 
J o 

dt[ ds(A lf + A 2f )u(t - s) Vl (s)A cl (0) + / ds(A lf + A 2f )u(t - s)v 2 (s)A d (t)} (A40) 

<ft[/ d Sl / (fs 2 / dS32(Ai / + A 2 /)7/(t-Si)[Gi2(si,S2)+G 2 l(52,Sl)] 

Jo Jo Jo 

xu(s 2 - s 3 )vi(s 3 )A d (0) (A41) 

+ d Sl ds 2 / ds 3 2(Ai/ + A 2 /)r ? (t-si)[Gi2(si,s 2 ) 
JO Jo Jo 

+G 2 i(s 2 ,si)]^(s 2 - s 3 )v 2 (s 3 )A d (i)]. (A42) 
To simplify the expressions, let us define 

= 2 / dsrf(t- s)u!(s), d 2 (t) = 2 / dsrf(t - s)u 2 (s) , (A43) 

JO JO 

ci(t) = dsi ds 2 / ds3r](t- si)[Gi 2 (si,s 2 ) + G 2 i(s 2 , si)Ms 2 - s 3 )«i(s 3 ), (A44) 
Jo Jo Jo 

c 2 (t) = ds 1 ds 2 / ds 3 r/(t - si)[Gi 2 (si,s 2 ) + G 2 i(s 2 ,si)]^(s 2 - s 3 )t; 2 (s3), (A45) 
Jo Jo Jo 

e iW = / dsi/(i — s)v 2 (s) = / — s)ui(t — s) = / dsu(s)ui(s) (A46) 

Jo Jo JO 

e 2 (i) = / dsi/(t — s)ui(s) = / dsu(t — s)u 2 (t — s) = / dsu(s)u 2 (s). (A47) 
Jo Jo Jo 

Finally, we have, 

J r {xif, x 2f ,y lf ,y 2f , t + dt\x w ,x 20 , y w , y 20 , 0) 

2 (.00 poo 

= N(t) I] / d/J** / d/? fej/ exp (— (Pi + ^ - - (3%)) 
x {1 - dt[i(V(x lf ,x 2f ) - V(y lf ,y 2f )) + i(A lf + A 2/ )(d 1 (t)(S i;1 + E i>2 ) 
+d 2 (t)(£i / + S 2/ )) + (Ai/ + A 2/ )(Ai,i + A i)2 )(e 2 (t) + 2ci(t)) 



+ (A 1/ + A 2/ ) 2 ( ei (t) + 2c 2 (t))]} 

2 [ dx\ f { l3lx) + dx\ f { ' ™> < dy ^ "V> < dy 2 f 



+ + ft(-^) 2 + S(-^) 2 + S("^) 2 ]}- (A48) 
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Hence 



J r (t + dt) -J r = -dtJ r [i^n 2 (x 2 lf + x 2 2f - y\ f - y 2 2f ) + {A lf + A 2/ ) x 

[i{d!(t)(Z hl + E i)2 ) + d 2 {t)(E lf + £ 2/ )) + (Ai,i + A ij2 )(e 2 (t) + 2ci(t)) 

+(A 1/ + A 2/ )(ei(t) + 2c 2 (t))]] (A49) 

ldtd 2 J r ldtd 2 J r ldtd 2 J r ldtd 2 J r 
+ 2~iBx{ f + 2~idxJ~ f + 2~dyY f + ^Tdyf/ { ' 

We can then get the evolution equation for the propagator J r . 

B 

i-^Jr(xif,x 2 f,yif,y2f, t\x 10 ,x 20 ,yi , y 20 , 0) 
B 

= i-^[Jr(xif,x 2 f,yif, y 2 f,t + dt\x w , x 2 o,y w ,y 20 , 0) - J r (x 1 f,x 2f , y lf , y 2 f,t\x 10 ,x 20 , y w ,y 2 o, 0)] 
1 B 2 B 2 B 2 B 2 1 

c ( U U a \ , r>2^ 2 , 2 2 2 \ 

= { " 2 { dxf f + dxff " ^ " dyj) * (Xlf X2f " Vlf " " 2/) 

+ (A lf + A 2/ )((di(t)(E M + E i>2 ) + d 2 (t)(Z lf + S 2/ ))) 

- i(A lf + A 2/ )(A M + A i)2 )(e 2 (t) + 2ci(i)) 

- i(Ai/ + A 2/ ) 2 (ei(t) + 2c 2 (t))} J r (xi/, x 2 /, yi/, y 2f , t\x w , x 20 ,y w ,y 20 , 0). (A51) 



4. Step four 

Now we have the explicit expression for J r . But we still need to deal with terms of the form 
like AuJ. To do so we can differentiate J with respect to Si/ and get 

d^ lf J = [Mt){A lf + A 2f )+ib 5 (t)(A lf 

-A 2/ ) - ib 3 (t)(A u + A 2l ) - ib 7 {t){A u - A 2i )]J. (A52) 

Similarly if we want A 2i J, we can differentiate J with respect to S 2 j and get 

5e 2/ J = [^i(t)(A 1/ + A 2/ )-z6 5 (t)(A 1J 

-A 2/ ) - ib 3 (t)(A u + A 2t ) + ib 7 {t)(A lt - A 2i )]J. (A53) 

The sum of these two equations gives 

(ds lf + a E2/ ) J = [MhWiAif + A 2/ ) - 2ib 3 (t)(A u + Am)] J. (A54) 
This can be written as 

(A u + A 2i ) J = t^—M^ + fc 2/ ) + 26i(t)(A 1/ + A 2/ )]J. (A55) 
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Similarly, we can differentiate with respect to A^ (or A 2 / )to get £iiJ (or £ 2 jJ). The sum of 
these two equations gives 

(dA lf +d A2f )J = 2[i& 2 (t)(Ei i + £ 2 i) + i6i(t)(£i / + E 2/ ) 

- oi 2 (t)(Ai i + A2i)-2oii(t)(Ai / + A 2 /)]J (A56) 

and 

(£i, + £ 2 *)J = ^[-i(9 Ai/ +9 A2/ ) + ^(a Ei/ +^ 2/ )-26 1 (*)(E 1/ + E 2/ ) 

-<[4fl U (t) + 2 °"yW + A 2 /)]^- (AST) 

Substituting in what we already have for (Eij + E 2 j) J and (Aij + A 2 j) J, and multiplying by 
and integrating over initial coordinates, we obtain 



(Ai / + A 2/ )di(t)(Ei i + E 2i )J 

_i /a ^a \> Q i2( f ) 



(A lf + A 2/ )d 1 (t)[^- 7 ^(a Al/ + d A2f ) + - + fc 2/ ) 



+ y ^ 2a n (t) Oi 2 (t)6i(t) 

MO (El/ + S2/) " * [ W 11 + 2/)] ' (A58) 



and 



(A i; + A 2/ )(e 2 (t) + 2c 1 (t))(A li + A 2i )J 

L 2b 3 (t)^ lf ' ' b 3 (t) 



= (A lf + A 2f )(e 2 (t) + 2c 1 (t))[^- 77T (9 El/ + ^ 2/ ) + AZ(A V + A 2/ )] J (A59) 



Hence we can write the evolution equation for the reduced density matrix as 

.8 . 1, ^ # # 9 2 . l o2 . 2 2 2 2,, 

^ = h 2 ( ^f + ^| " " W 2 ] + 2 {Xl + X2 " Vl " " 2)K 

+«2 2 (t)(A 1/ + A 2/ )(E 1/ + E 2/ ) Pr 
-iAi(t)(Ai/ + A 2f )(d Alf + d Aa/ )pr 
-iA 2 (t)(Ai / + A 2/ ) 2 p r 

+A 3 (*)(A 1/ + A 2/ )(d Sl/ + a S2/ )/> r (A60) 



where 



5 - 5 5 ■ 5 _ 1 d 5 A 



and 



<50 2 (t) = d 2 (t)- dl (t)M|, ^(i)^^, (A62) 
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Mt) - m\^§ + ] + (*(*) + *>(*)) + (*(*) + 2* (*))£$ (A63) 

_ di(t)a 12 (t) e 2 (t) + 2 Cl (t) 

^ = + 263(«) • (A64) 

This immediately leads to the general master equation 



5. Coefficients of the Master Equation 

The determination of the coefficients is reasonably standard, so we only provide the explicit 
forms of those time-dependent functions that will be used later on. As shown in [ljj, the functions 
5Q 2 (t),T(t), A(i), £(t) can be constructed in terms of the elementary functions Ui(s),i = 1, 2, which 
satisfy the following homogeneous integro-differential equation: 

f(s) + n 2 f(s) + -i J* d\ V (s - X)f(X) = (A65) 

with the boundary conditions: 

Ul (s = 0) = 1 , m(s = t) = 0, (A66) 

and 

n 2 (s = 0) = , u 2 {s = t) = l. (A67) 
Here rj(t — s) is the dissipation kernel given by 

POO 

rj(s) = — / du}I(uj) sm(u!s), (A68) 
Jo 

and I(uj) is the spectral density of the environment. Note that the numerical factor 4 before the 
integral in this equation is different from that in . This is the main difference due to the presence 
of two harmonic oscillators. Although the two harmonic oscillators are not coupled directly, they 
are connected by the common reservoir, hence they affect each other dynamically. 
Let G\(s,t) be the Green function obeying the following equation: 
;2 



^Gi( S ,r) + n 2 Gt(s,T) + ^ j^drvis - t)Gx(s,t) = 6(s - r), 



(A69) 



with initial conditions: 



Gi(s = 0,t) = 0, ^-Gi(a,r)| a=0 = 0. (A70) 
as 
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The Green function G2(s,t) is defined analogously. The coefficients can then be written as 

^ (t , = ^ (t -.)(^)-2^), (ATI) 

m = ±- fdsvit-s)^-, (A72) 
M Jo ui(t) 

A{t) = 2^£ dAGl( *' AM *" A) 

~w J* ds Jl dT Jo dXr]{t ~ s)Gl(t ' A)G2(s ' rMr ~ A) ' (A73) 



and 



£(t) = - / d\G[(t, \)v(t - A) 

Jo dS J dT Jo dXV ^ ~ ^i^' A ) G 2(s,t)z/(t - A), (A74) 

where defined as 

/•+00 J 

i/(s) = / dooI(uj) coth.(-Huj(3) cos(us), (A75) 

Jo 2 

is the noise kernel of the environment. Here a "prime" denotes taking the derivative with respect 
to the first variable of G\(s,t). 

APPENDIX B: EXPLICIT EXPRESSIONS FOR Pij 

We find that the matrix Gij is the same for all the p^. Thus, we can write Gij = G. The matrix 
elements for the matrix G are given by 



Gn 


= G22 


G33 


= G44 




= G21 


G34 


= G43 


G13 


= G14 



z&4 ibs 1 

ft22 + T + T + W 

ib^ ibs 1 

-(2a 22 + ibi - ib & ), 
^(2a 2 2 - ibi + ibs), 

G14 = G23 = G24 = G31 = G32 = G41 = G42 = — 0>22- (B5) 



(Bl) 
(B2) 
(B3) 
(B4) 



Then the determinant of G can be explicitly computed, 

1 Q22 b\ 
165 8 25 6 25 A ' 4(5 4 ' 5 2 ' 



det G = b A b s + —Ts + TTIfi + TIT + 77T + ~ ra" ■ ( B6 ) 
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Moreover, the matrix elements of the inverse matrix G 1 are 



ihh 2 xj_4. __J^ _ ___t _ __^L _l 3_ _ ifl 22 & 8 _§L rR7 x 
"2 4 8 + 8<5 6 4<5 4 85 4 85 4 A5 2 5 2 + 4<5 2j ' 1 ' 



G 33 ~ G AA - ditG^2^ 8 + ° 22 ^ 



+ 2 4&8 + 8F + lF + 8^ + 8^ + 4^ + "T 2 "" + W>' (B8) 



Gr^G^ 1 = ^a^\b, + a 22 bl 



Jhh 2 _ ° 22 _ i&4 a- i&8 6 4 , ^22 &8 &1 x /r>Q\ 

2 &4 ° 8 4<5 4 8<5 4 + 8<5 4 4<5 2 + 5 2 + 45 2j ' 1 9J 

C34 = G43 1 = ^(-l & 2 6 8 + a 22 6 2 

Jhh 2 ° 22 , i6 4 i&s ^22 &8 , bj 

+ 2 Ms "4^ + 8^"8^"4^" "T 2 "" + AP h (BW) 

G 13 = G U G 23 = G 2A = G 31 = G 32 = G A1 = G A2 = ^Q^ 22 ^ + ^a)' ( BU ) 

For the case of p\\: 

f . n x AT A r (glO ~ L 0) 2 + (^20 ~ A )) 2 + (2/10 ~ Lq) 2 + (j/ 2 Q ~ £ ) 2 1 

pn(t = UJ = iv exp [ 



2<5 2 J 
x exp[«P (a;io + X20 -yio - 2/20)], (B12) 



then the matrix elements for F are, 









iP - a\ 2 xi + 


ib 2 x\ 
2 


z& 3 xi i& 6 xi 
2 2 






+ 


ib 2 yi 


«&3J/1 

2 


ihyi ihyi 
2 2 




= iP 




ib 2 xi 
a\2X\ H — 


2 


i&6Xi ib-jXi 
2 ' 2 






+ 


ib 2 yi 
a\2Vi 2 


ibiyi 
2 


. ihyi , ^72/1 
+ 2 ' 2 




-iP 


+ 


ife 2 xi 

«12Xl + 


2 


2 2 








0121/1 2 


^32/1 
2 


ihyi ihyi 
2 2 


Ftl = 


-iP 


+ 


i6 2 xi 

012^1 + 


+ 2 


i&6Xi ihx\ 
2 2 








^2^1 

«i2yi 2 


^32/1 
2 


2 2 



^67X1 *t>2^2 i&3^ 2 ihx 2 ibjx 2 
- — - aM + ~2 2 2" + — 

*&22/2 ^32/2 , i£>62/2 , «&72/2 , ^0 

+ a 12 y 2 -— — + — + — + - 

ib 2 x 2 ib 3 x 2 ib 6 x 2 ib 7 x 2 
■ - a\ 2 x 2 H 1 — 



ib 2 y 2 ihy 2 ib 6 y 2 ib 7 y 2 L 
+ a 12 y 2 -— - - 2" + ^ 

ib 2 x 2 ib 3 x 2 ib 6 x 2 ib 7 x 2 
■ + ai 2 x 2 H 1 — 



ib 2 y 2 , ihV2 , ihV2 ihV2 , L 

- a 12 y 2 ~ — + — + — 2 2~ + ^ 

ib 2 x 2 ib 3 x 2 ib 6 x 2 ib 7 x 2 
+ a l2 x 2 + — + — + — + — 

ib 2 y 2 , ihy 2 ihV2 , ihV2 , 

- a 12 y 2 - — + — _ + ^ + _, 
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where f£ = (ift , F 2 , , ift , ift ) and 

Q Z Z Q 

cn = — anXj + 2 &lX i + 2^ 5Xl ~~ 2ana;iX2 + ib\X\x 2 - ib 5 x\x 2 

Z Z 

- aux 2 + -hxl + —65^1 + 2onXiyi + 2a n x 2 yi 

~ aiiVi ~ ~ 2 &5y i + 2a n x iy2 + ib\xiy 2 + 2a u x 2 y 2 + ib\x 2 y 2 

i i 2L 2 

- ZaiiVm + ihvm - a\\y\ + ^hy 2 - -hyi - (B13) 

For the case of p\ 2 : 

n x AT 4 r (SIP ~ L o) 2 + (^20 ~ Lp) 2 + (gig - Lp) 2 + (y 20 + Xp) 2 1 
Pi2{t = U) = iV exp [ ^2 J 

x exp [iP (x 10 + x 20 - 2/10 + 2/20)] (B14) 

Fi2=Fn, *i2 = *n, *i2 = *ii, F? 2 = F? 1 + 2iP -2^, c 12 = c n . (B15) 

For the case of P13: 

n x W 4 r ( g io ~ L o) 2 + (^20 - Lp) 2 + (yig + Lp) 2 + (3/20 - X ) 2 1 
Pi3(t = UJ = iV exp [ ^2 J 

x exp [iPp(xip + x 20 + yio - y 20 )\ (B16) 

^13=^11, *i3 = *ii, ^3 = ^1 + 2^0-2^, i? 3 = *n, c 13 =c n . (B17) 

For the case of p\±: 

,. n , AT 4 r (&10 - ^o) 2 + (^20 - ^o) 2 + (yio + U) 2 + (y2o + L ) 2 , 
Pl4 (t = 0) = iV exp[ ^ ] 

x exp [iP (x w + x 2 p + yw + 2/20)] (B18) 

Fh = Fl 1 , F& = F? 1 , F? 4 = F? l + 2iP -2^-, F? 4 = + 2iP - 2^, c 14 = c n (B19) 
Similarly, one can work out the cases for p 2 i, p%i and pa (i = 1, 2, 3,4). 
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